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11 Abstract 
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12 Recent studies, both based on remote sensed data and coupled models, showed a re- 

13 duction of biological productivity due to vigorous horizontal mixing in upwelling systems. 
I— I 14 In order to better understand this phenomenon, we have considered a system of oceanic 
^ 15 flow in the Benguela area coupled with a simple biogeochemical model of Nutrient-Phyto- 
^ 16 Zooplankton (NPZ) type. For the flow three different surface velocity fields are considered: 
^ 17 one derived from satellite altimetry data, and the other two from a regional numerical 
^ 18 model at two different spatial resolutions. We computed horizontal particle dispersion 
^ 19 in terms of Lyapunov Exponents, and analyzed their correlations with phytoplankton 

> 20 concentrations. Our modelling approach confirms that in the south Benguela, there is 

21 a reduction of biological activity when stirring is increased. Two-dimensional offshore 

22 advection seems to be the dominant process involved. In the northern area, other factors 

23 not taken into account in our simulation are influencing the ecosystem. We provide expla- 

24 nations for these results in the context of studies performed in other Eastern Boundary 

25 upwelling areas. 
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26 1. Introduction 



Marine ecosystems of tlie Eastern Boundary Upwelling zones are well known for their 
major contribution to the world ocean productivity. They are characterized by wind- 
driven upwelling along the coast of cold nutrient-rich waters that supports elevated plank- 



30 ton and pelagic fish production (Mackas et al. , 2006). Variability is introduced by strong 



advection along the shore, physical forcings by local and large scales winds, and high sub- 
and mesoscale activities over the continental shelf and beyond, linking the coastal domain 
with the open ocean. 

The Benguela Upwelling System (BUS) is one of the four major Eastern Boundary 
Systems (EBUS) of the world. The coastal area of the Benguela ecosystem extends from 
southern Angola (around 17°S) along the west coast of Namibia and South Africa (36°S). 
It is surrounded by two warm temperate boundary currents, the Angola Current in the 
north, and the Agulhas Current in the south. The BUS can itself be subdivided into two 



subdomains by the powerful Luderitz upwelling cell (Hutchings et al. , 2009). Most of the 



biogeochemical activity occurs within the upwelling front and the coast, although it can 
be extended further offshore toward the open ocean by the numerous filamental structures 



developing offshore (Monteiro, 2009). In the BUS, as in the other major upwelling areas. 



a high mesoscale activity due to eddies and filaments is observed and impacts strongly on 



marine planktonic ecosystem over the shelf and beyond (Brink and Cowles, 1991 Martin 



2003 Sandulescu et al , 2008 Rossi et al. , 2009). 



The purpose of this study is to analyze the impact of the horizontal stirring on the phy- 



toplankton dynamics in the BUS. Recently, Rossi et al. (2008, 2009), using satellite data 



of the ocean surface, suggested that mesoscale activity has a negative effect on chlorophyll 
standing stocks in the EBUS. This was obtained by correlating remote sensed chlorophyll 
data with a Lagrangian measurement of lateral stirring in the surface ocean (see Methods 
section below). This result was unexpected since mesoscale transport, particularly due 
to eddies, has been related to higher planktonic production and stocks in the open ocean 



2 



53 (McGillicuddy et al. , 2007) as well as off a major EBUS ( Correa- Ramirez et al. , 2007). 



A more recent and thorough study performed by Gruber et al. (2011) in the California 



and the Canary current systems detailed the initial results from Rossi et al. (2008, 2009). 



Based on satellite derived estimates of net Primary Production, of upwelling strength and 
of Eddy Kinetic Energy (EKE) as a measure the intensity of mesoscale activity, they con- 
firmed the suppressive effect of mesoscale structures on biological production in upwelling 
areas. The mechanism behind this observation was investigated using 3D eddy resolving 
coupled models. The eddies tend to export offshore and downward a certain pool of nu- 
trients not being effectively used by the biology in the coastal areas. This process they 
called "nutrients leakage" is also having a negative feedback effect by diminishing the 
nutrients available in the deep waters being re-upwelled continuously. 

In our work, we focused on the Benguela area, being the most contrasting area of 
all EBUS in term of mixing intensity. Although mechanisms involved occur in the 3D 
space, the initial observation of this suppressive effect was based only on two-dimensional 



(2D) datasets (Rossi et al. , 2008). Here we use 2D numerical analysis in a simple semi- 
realistic framework to test the effect of horizontal advection versus biological dynamics. 
Meanwhile, since vertical dimension is crucial in upwelling areas, it was introduced in 
our model in a simplified way by considering a source term with an intensity and spatial 
distribution corresponding to the upwelling characteristics. Indeed other theoretical stud- 
ies in idealized 2D settings display also negative correlation between mixing and biomass 



73 (Tel et al. 2005; MacKiver and Neufeld 2009). Contrarily to EKE which is an Eulerian 



diagnostic tool, we used here a Lagrangian measurement of mesoscale intensity. It has 
been demonstrated as a powerful tool to study patchy chlorophyll distributions due to 



76 dynamical structures at mesoscale, such as upwelling filaments (Calil and Richards, 2010). 



Different velocity fields were considered, one obtained from satellite and others from nu- 
merical simulations. The robustness of our results with respect to spatial resolution is 
tested by using two numerical velocity datasets at different resolution. Our results are 
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80 compared with real distributions of chlorophyll (a metric for phytoplankton) obtained 

81 from SeaWiFS satellite. 

82 This paper is organized as follows. In Sec. |2] we describe the different data sets 

83 for our analysis. Sec. [3] contains the methodology, including the Finite-Size Lyapunov 

84 exponents, and the plankton numerical model. Then, in Sec. |4]our results are presented 

85 and discussed in the context of existing bibliography. Finally in Sec. [5} we summed-up 

86 our main findings. 



87 2. Satellite and simulated data. 

88 A total of three sources of two-dimensional velocity data sets in the surface of the 

89 Benguela area were used: two were obtained from the numerical model ROMS (Regional 

90 Ocean Model System), and the other one from a combined satellite product. ROMS is 

91 a free surface, hydrostatic, primitive equation model, and the run used here was eddy 



resolving but climatologically forced (Gutknecht et al. 2011). At each grid point, linear 



93 horizontal resolution is the same in both the longitudinal, 0, and latitudinal, 9, directions, 

94 which leads to angular resolutions A0 = Aq and = A0cos0. The numerical model 

95 was run onto 2 different grids: a coarse one at Aq = 1/4°, and a finer one at Aq = 1/12° 

96 of spatial resolution. In the following we label the data set from the coarser resolution 

97 as ROMSl/4, and the finer one as ROMSl/12. In both of them, vertical resolution is 

98 variable with 30 layers in total. Only data from the upper layer were used. The third 

99 set of velocity data are surface currents computed from a combination of wind-driven 

100 Ekman currents, at 15 m depth, derived from Quickscat wind estimates, and geostrophic 

101 currents calculated using time variable Sea Surface Heights (SSH) obtained from satellite 



102 (Sudre and Morrow, 2008). These SSH were calculated from mapped altimetric sea level 

103 anomalies combined with a mean dynamic topography. This velocity field, labeled as 

104 Satellitel/4, covers a period from June 2002 to June 2005 with a spatial resolution of 

105 Ao = 1/4° in both longitudinal and latitudinal directions. 
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106 To validate simulated biological fields we used a tliree-year-long time series, from Jan- 

107 uary 2002 to January 2005, of ocean color data. Phytoplankton pigment concentration 

108 (chlorophyll-a) are obtained from monthly SeaWiFS (Sea viewing Wide Field-of-view Sen- 

109 sor) products, generated by the NASA Goddard Earth Science (GES) /Distributed Active 
no Archive Center (DAAC). Gridded global data were used with a resolution of approxi- 
111 mately 9 by 9 km. 



112 3. Methodology. 



113 3.1. Finite Size Lyapunov Exponents (FSLEs). 



FSLEs ( Artale et al. 1997 Aurell et al. 1997; Boffetta et al. 2001 ) provides a measure 



of dispersion, and thus of stirring and mixing, as a function of the spatial resolution, 
serving to isolate the different regimes corresponding to different length scales of the 
oceanic flows, as well as identifying the Lagrangian Coherent Structures (LCSs) present 
in the data. FSLE are computed from r, the time required for two particles of fluid (one 
of them placed at x) to separate from an initial (at time t) distance of to a final distance 
of (5/, as 

A(x,t,5o,5/) = -log^. (1) 
r do 

It is natural to choose the initial points x on the nodes of a grid with lattice spacing 
coincident with the initial separation of fluid particles 5q. Then, values of A are obtained 
in a grid with lattice separation 5o- In this work we take always the resolution of the FSLE 
field, equal to the resolution of the velocity field, Aq. Other choices of parameter 
are possible and 5q can take any value, even much smaller than the resolution of the 



velocity field (Hernandez-Carrasco et al. , 2011a). This opens many possibilities that will 



not be explored in this work, since we focus here in the primary production, and, in 
some instances, the influence of the data resolution, not on the resolution of the FSLEs 
computation. 
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123 The field of FSLEs thus depends on the choice of two length scales: the initial, 



ip,nd the final 5f separations. As in previous works (d'Ovidio et al. , 2004, 2009 Rossi 



et al. , 2008 Hernandez- Carrasco et al. , 2011a) we will focus on transport processes at 



mesoscale, so that 5f is taken as about 110 km, which is the order of the size of mesoscale 
eddies at mid latitudes. To compute A we need to know the trajectories of the particles 
which gives Lagrangian character to this quantity. The equations of motion that describe 
the horizontal evolution of particle trajectories in longitudinal and latitudinal spherical 
coordinates, x = (0, A), are: 



dcf) 
~dt 

de 

~dt 



Rcos9 
R ' 



(2) 
(3) 



where u and v represent the eastwards and northwards components of the surface velocity 
field, and R is the radius of the Earth (6400 km). 

The ridges of the FSLE field can be used to define the Lagrangian Coherent Struc- 



tures (LCSs) (Haller and Yuan, 2000 d'Ovidio et al. 2004 2009; Tew Kai et al. 2009 



Hernandez-Carrasco et al. , 2011a), useful to characterize the flow from the Lagrangian 



136 point of view (Joseph and Legras, 2002 Koh and Legras, 2002). In fact, since we are only 



interested in the ridges with large values of FSLE, the ones which significantly affect mix- 
ing, LCSs can be obtained as the regions with high values of FSLE, which have a line-like 
shape. We will compute FSLEs integrating backwards-in-time the particle trajectories, 
since attracting LCSs associated to this (the unstable manifolds) have a direct physical in- 
terpretation (Joseph and Legras 2002; d'Ovidio et al. , 2004 2009[ ). Tracers (chlorophyll, 
temperature, ...) spread along the attracting LCSs, thus creating their typical filamental 



structure (Lehan et al. , 2007 Calil and Richards 2010). 



144 3.2. The Biological model 



145 - The plankton model is similar to the one used in previous studies by Oschlies and 



146 Gargon (1998, 1999) and Sandulescu et al. (2007, 2008). It describes the interaction of 
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147 
148 
149 



a three-level trophic chain in the mixed layer of the ocean, including, phytoplankton P, 
zoo-plankton Z and dissolved inorganic nutrient A^, whose concentrations evolve in time 
according to the following equations: 



= = - /3— + /i^ (1 - 7) + /ipP + /i.^' , (4) 
dt + N \ a + TjP^ J 

dP ^ N ariP^ 

dt kn + N a + rjP^ 

150 where the dynamics of the nutrients, Eq. Q, is determined by nutrient supply due to 

151 the vertical mixing $jv, its uptake by phytoplankton (2"'^ term) and its recycling by 

152 bacteria from sinking particles term). Vertical mixing which brings nutrients from 

153 lower layers into the mixed surface layer of the ocean is parameterized in the model (see 

154 below), since the hydrodynamical part considers only horizontal 2D transport. Terms in 

155 Eq. ([5]) stand for phytoplankton growth by consuming A^, the grazing by zooplankton, 

156 and its natural mortality. The last equation, Eq. (|6]), represents zooplankton growth by 

157 consuming phytoplankton minus its quadratic mortality. 

A crucial part of this model comes in the vertical mixing, <l>jv, since it mimics the 
upwelling. Assuming constant nutrient concentration Aq below the mixed layer, this 
term reads: 

$^ = 5(x,t)(A^o-iV), (7) 

158 where the temporally and spatially dependent (on the two dimension location x) function 

159 5* determines the strength and the horizontal spatial distribution of vertical mixing in 

160 the model, thus specifying the upwelling characteristics. Thus, the vertical dynamics is 

161 introduced in our two-dimensional model via this function S. Upwelling intensity along 

162 the coast is characterized by a number of cells of enhanced vertical ekman driven transport 



163 that are associated with similar fluctuations of the alongshore wind (Demarcq et al. , 2003 



Veitch et al. , 2009). Following these results, we use a function S which is different from 
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zero in a strip 0.5° wide from the coast. Its spatial dependence along the coast is plotted 
in Fig. [T] For the temporal dependence, S alternates between the two configurations 
displayed in Fig. [T]one for winter and another for summer. Six separate upwelling cells can 
be discerned in the figure, with peaks at approximately 33°S, 31°S, 27.5°S, 24.5°S, 21.5°S, 
17.5°S, which are known with the following names: Peninsula, Columbine+Namaqua, 
Luderitz, Walvis Bay, Namibia and Cunene, respectively. Luderitz being the strongest. 

The dynamical system given by Eqs. ( ipp ), for values of S in the range shown on 
Fig. [T| evolves towards equilibrium for A^, P and Z. But S is not fixed and its spatial 
dependence introduces a coupling with the hydrodynamics. The transient time to reach 



equilibrium is typically 60 days with the initial concentrations used (see Sec. 3.3). The 



parameters are set following a study by Pasquero et al. (2004) and are listed in Table [T] 



parameter 


value 






0.66 day~^ 




T) 


1.0 (mmol N m~^)~^ 


day ^ 


7 


0.75 




a 


2.0 day-i 




kN 


0.5 mmol N m~^ 






0.2 




jJLp 


0.03 day-i 






0.2 (mmol N m-^)-^ 


day~^ 


No 


8.0 mmol N m~^ 





Table 1: List of parameters used in the biological model. 



176 3.3. Coupling hydrodynamical and biological model in Benguela. 

177 The evolution of the concentrations within a fiow is determined by the coupling be- 
ns tween the hydrodynamical and biological models, and it is performed by the advection- 



8 



0.5 1 — \ — ' — \ — ' — \ — I — \ — ' — \ — ' — \ — ' — \ — ' — \ — ' — \ — I — [ 




I I I I I I I I I I I I I I I I I I I 

34S 32S 308 28S 26S 24S 22S 20S 18S 16S 
Latitude -0- (degrees) 



Figure 1: Shape and values of the strength (S) of the upwelhng cells used in the simulations for winter 



and summer seasons (following Veitch et al. (2009)). 



179 reaction-diffusion system. Thus, the complete model is given by the following system of 

180 partial differential equations: 



ON 

~dt 
dP 

~dt 

dz^ 

~dt 



+ vVAT = Fn + DV^N, 
vVP = Fp + DV^P, 
vVZ = Fz + DV^Z. 



(9) 
(10) 



181 The biological model is the one described before by the functions F/v, Fp and Fz- 

182 Horizontal advection is the 2D velocity v, which is obtained from satellite data or from 
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the ROMS model. We add also an eddy diffusion term, via the operator, acting on 
A^, P, and Z to incorporate the small-scale turbulence, which is not explicitly taken into 
account by the velocity fields used. 



The eddy diffusion coefficient, is given by Okubo's formula (Okubo, 1971), -D(0 = 
2.055 * 10~^ Z^'^^, where / is the value of the resolution, in meters, corresponding to the 
angular resolution / = Aq. The formula gives the values Z)=26.73 vn? j s for Satellitel/4 
and ROMS 1/4, and D=7A m^s for ROMS 1/1 2. 

The coupled system Eqs. ( SfOpIO ) is solved numerically by the semi-Lagrangian algo- 



191 rithm described in Sandulescu et al. (2007), combining Eulerian and Lagrangian schemes. 



The initial concentrations of the tracers were taken from Kone et al. (2005) and they 
are A^o = 1 TnmolNm~^ , Pq = 0.1 mmolNm~^ , and Zq = 0.06 mmolN m~'^ . The in- 
flow conditions at the boundaries are specified in the following way: into the eastern, 
western, and southern parts of the computation domain fluid parcels enter with very 
poor biomasses concentration: Nl = O.OlA^^o mmolNm~^ , Pl = O.OIPq mmolNm~^, and 
Zl = O.OIZq mmolNm"^ . Across the northern boundary, fluid parcels enter with higher 
concentrations Nh = 5 mmolNm~^ , Ph = 0.1 mmolNm~^, and Zh = 0.06 mmolNm~^ 



199 according with the values given by CARS for the Benguela system (Condie and Dunn 



2006). The integration time step is (it = 6 hours. 



201 4. Results and discussion. 

202 In this section we first compute the FSLEs on the velocity fields to quantify the 

203 horizontal stirring activity over the area. Then we analyze the results of the coupled 

204 biological-hydrodynamic model. Finally we investigate the relation between horizontal 

205 stirring activity and biological productivity. 

206 Horizontal activity 

207 We have computed the FSLE with a initial separation of particles equal to the spatial 

208 resolution of each velocity fields (5o= 1/4° for Satellitel/4 and ROMSl/4, and 6^= 1/12° 
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209 for ROMSl/12). As already mentioned, the final distance is always chosen to focus on 

210 transport processes by mesoscale structures at mid latitudes, 5f= 1°. The areas of more 

211 intense horizontal mixing can be identified using time averages of the backward FSLEs 



212 (d'Ovidio et al. , 2004). Figure |2j allows an easy characterization of sub-regions with dif- 

213 ferent horizontal mixing activity in the Benguela system. Areas of large average values of 

214 FSLEs are identified as exhibiting an intense horizontal stirring or mesoscale activity. We 



215 confirm the results of Rossi et al. (2009) by using different velocity data sets. Although 

216 there are visible differences in the detailed patterns, good agreement between all datasets 

217 is shown when computing the spatial correlation: for instance, correlation coefficient E? 

218 between FSLEs map from Satellitel/4 and from ROMSl/4 is 0.81. Correlation coeffi- 

219 cients between Satellitel/4 and ROMSl/12 on one hand, and between ROMSl/4 and 

220 ROMSl/12 on the other hand, are lower (0.61 and 0.77 respectively) since the FSLE were 

221 computed on a different resolution. More details on the effect on the grid resolution when 



222 computing FSLEs can be found in Hernandez- Carrasco et al. (2011a). For all data sets 

223 high mixing values are observed in the southern region, while the northern area displays 

224 significantly lower values. Note that the separation is well marked for Satellitel/4 where 

225 the line between the two areas is around 27°. In the case of the ROMS data sets, the 

226 mixing activity is more homogeneously distributed, although the north-south gradient 

227 is still present. We associate this difference with the injection of strong and numerous 

228 Agulhas rings into the south of the area from the Agulhas retroflection. 

229 The latitudinal behavior of mixing along the coastal upwelling can be seen in Fig. |3| 

230 This was performed by computing the longitudinal averages of the plots in Fig. [2] for 

231 two coastally oriented strips, of 3° and 6° width, respectively. It is clear that horizontal 

232 mixing decreases as latitude decreases. Note that there are differences in the mixing 

233 values (FSLEs) depending on the type of data, their resolution and the grid size of FSLE 

234 computation. In general, considering velocities with the same resolution, the lower values 

235 correspond to Satellitel/4 as compared to ROMSl/4- On average, values of mixing from 
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ROMS 1/4 are larger than those from ROMS 1/1 2, whereas we would expect the opposite 
considering the higher resolution of the latter simulation favouring small scales processes. 
However a caveat here is that FSLE were not computed on the same resolution, so there 
are not directly comparable. Note also that a low-mixing region is observed from 28° to 
30° S on all calculations. It seems to indicate that the ROMS model is representing pretty 



241 well the spatial variability of the mixing. As proposed in a recent study by Titaud et al. 



242 (2011), these preliminary results indicate that FSLEs could be used as a diagnostic to 



validate eddy-resolving oceanic models. 

In Fig. |3] (bottom) we see that, for Satellitel/4, the values of FSLEs decay from 
0.18 days~^ in the southern to 0.03 days^^ in the northern area, with similar decays 
for ROMS 1/4. Specifically the North-South difference for Satellitel/4, ROMS 1/4 and 
ROMSl/12 are of the order of 0.15 days~^ , 0.15 days~^ and 0.08 days^^, respectively, 
confirming a lower latitudinal gradient for the case of ROMSl/12. These values do not 
change much when it is averaged over the 3 degrees stripe offshore (Fig. |3| top), although 
in this case relative maxima and minima appear, probably in relation with the complex 
and variable shelf circulation. 

The mixing behavior can be also assessed by looking at a proxy of the intensity of 



mesoscale activity, the Eddy Kinetic Energy (EKE), as done in Gruber et al. (2011 ). Fig. E 



shows that there are regions, as in the FSLE case, with distinct dynamical characteristics. 
Larger values appear in the south and smaller in the north. This distribution is in good 
agreement with the one deducted from the FSLEs (Fig. |2]). Some simple spatial correla- 
tion (not shown) indicate that EKE and FSLE patterns are well correlated when using a 
non-linear fitting (power law). For instance, EKE and FSLE computed on the velocity field 
from SatelUtel/4 exhibit a of 0.86 for the non-linear fitting: FSLE = 0.009- EKE^-^^ . 



260 It is in agreement with the initial results from Waugh et al. (2006); Waugh and Abraham 



261 (2008), for a related dispersion measurement, and confirmed the thorough investigation 



262 of the relationship between EKE and FSLE by Hernandez- Carrasco et al. (2011b). 
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263 In the following sections, we study the effect of this variable surface mixing activity 

264 on the plankton dynamics. 



265 4-^- Plankton dynamics in the Benguela upwelling system. 

266 Evolution of A^, P and Z over space and time is obtained by integrating the systems 

267 described by Eqs. 8][9p!0 The biological model is coupled to the velocity field after the 

268 transient time needed to reach stability (60 days). In Figjsjwe show some snapshots of 

269 phytoplankton concentrations for the three velocity fields at different times. Since both 

270 ROMS simulation were climatologically forced runs, the dates do not correspond to a 

271 specific year, whereas we used the actual date for Satellitel/4- The most relevant feature is 

272 the larger value of concentrations near the coast due to the injection of nutrients following 

273 Fig. [T| Obviously the spatial distribution of P is dominated by the submeso- and meso- 

274 scale structures such as filaments and eddies. This is specially noticeable in the south, 

275 due to the presence of several Agulhas rings, cyclonic eddies and filaments. Differences 

276 are however observed for the three data sets. In particular, it seems that for Satellitel/4 

277 and ROMSl/12 the concentrations extend farther offshore than for ROMSl/4- 



278 Several studies (Lehan et al. 2007 d'Ovidio et al. , 2009 Calil and Richards 2010) 

279 have shown that chlorophyll distributions in the marine surface are linked to the local 

280 maxima or ridges of the FSLEs. This also occurs in our numerical setting, as it is visually 

281 shown in Fig. [6] We superimpose contours of high values of FSLE (locating the LCS) 

282 on top of phytoplankton concentrations for ROMSl/12 (every 8 days during a 32 days 

283 period). In some regions P concentrations are constrained and stirred by lines of FSLE. 

284 For instance, the edges of the cyclonic and anti-cyclonic eddies centered at 6 °E, 32 °S, 

285 and 28 °S in Fig. [6] on June 11 exhibit large values of phytoplankton concentration. This 

286 reflects the fact that tracers, even active such as chlorophyll, still disperse along these 

287 LCSs. 

288 In order to reveal regions of more intense biological activity, we have computed the 

289 temporal average of simulated P. The results, plotted in FigjT] a), b), c), show that 
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290 coastal regions with high P extend approximately, depending on latitude, between half 

291 a degree and two degrees offshore. It is comparable with the pattern obtained from the 

292 satellite-derived chlorophyll data (FigjTjd)). The spatial correlation of averaged simulated 

293 chlorophyll with satellite is as follows: = 0.85 for Satellitel/4 versus SeaWIFS; B? = 

294 0.89 for ROMSl/4 versus SeaWIFS and = 0.85 for ROMS 1/1 2 versus SeaWIFS. 

295 Despite the very simple setting of our models, the phytoplankton development over the 

296 Benguela shelf is well simulated by the upwelling parameterization chosen. Note however 

297 that our simulated chlorophyll values are about ~ 3-4 times lower than satellite data, 

298 as shown by the colorbar scale. Of course several factors, both biological and physical, 

299 are not taken into account in this simple setting that might explain this offset. Another 

300 possible explanation is the low reliability of the ocean color in very coastal waters optically 

301 complex. 

302 We now examine the latitudinal distribution of P. The top row in FigjS] displays the 

303 outputs of the numerical simulations that were averaged over a coastal strip of 3° (left) 

304 and 6° (right) width. The bottom row is the same but from the satellite chlorophyll data. 

305 First of all, phytoplankton biomass has a general tendency to decrease with latitude, an 

306 opposite tendency to the ones exhibited by mixing (from FSLEs and EKE) for the three 

307 data sets. P values are higher in the northern than in the southern area of Benguela. A 

308 common feature is the minimum located just below the Luderitz upwelling cell (28°S), 

309 maybe related to the presence of a physical boundary, already studied and named the 



310 LUCORC barrier by [Shannon et al.| ( |2006D and [Lett et al.| ( |2007D . Note that on Fig. [3j 

311 (upper plot), the same latitude was marked by a local maximum of mixing that might 

312 be responsible for this barrier. Though not so evident, the same latitudinal tendency is 

313 observed for the SeaWIFS data plotted in Fig. ^) and d). Correlation of zonal average of 

314 simulated chlorophyll versus satellite data does not give striking results when considering 

315 the whole area {R^ ranging from 0.1 to 0.5). However, when considering each subsystem 

316 independently, high correlation coefficients are found for the south Benguela (i?^ around 
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317 0.75), but not for the north. It clearly indicates that our simple modelling approach is able 

318 tostimulate well the spatial patterns of chlorophyll in the south Benguela, but not properly 

319 in the northern part. The 2D vigorous mixing in the south and its associated intense off- 

320 shore export are sufficient to explain reasonable latitudinal patterns of P. The numerous 

321 eddies released from the Agulhas system, moving offshore in the south Benguela, might 

322 limit the large development of P by exporting unused nutrients and young phytoplankton 



323 communities toward the open ocean, as stated by Gruber et al. (2011). It also suggests 

324 that the negative effect seems to be mainly driven by 2D advection toward the open 

325 ocean. In the north, other factors seem to play an important role. Among many others, 

326 the 3D flow, the shelf width, the rivers and aeolian inputs, the remineralisation pattern, 

327 the presence of particular biogeochemical functioning,... etc. have been disregarded from 

328 this study, whereas they seem to impact widely plankton dynamics in the north. 

329 To address the question of the negative effect of horizontal stirring on phytoplankton 

330 concentration in a more quantitative way, we have examined the correlation between 

331 these two quantities. We have plotted spatial averages over each subregion (North and 

332 South) of every weekly map of FSLE versus the same average of the corresponding weekly 

333 map of P, for each week during three years in the case of Satellite and for one year for 

334 the case of ROMS (Figjoj). For all cases, a negative correlation between FSLEs and 

335 chlorophyll emerges. Thus, the higher the surface stirring/mixing, the lower the biomass 

336 concentration. The correlation coefficient is quite similar for all the plots (i?^=0.80 to 

337 0.84), and the slopes have the following values: -1 for Satellitel/4, -0.65 for ROMSl/4 



338 and -1.5 for ROMSl/12. Note that, similarly to the results of Rossi et al. (2008, 2009) 



339 and Gruber et al. (2011), the negative slope is larger but less robust when considering the 

340 whole area rather than within every subregion. The suppressive effect of mixing might be 



341 dominant only when mixing is intense, as in the south Benguela. Moreover, Gruber et al. 



342 (2011 ) stated that the reduction of biomass due to eddies may extend beyond the regions 

343 of the most intense mesoscale activity, not considered here. In fact in our simulations, we 
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344 observe than averaging FSLEs over a 3° or a 6° coastal band returns quite comparable 

345 absolute values, attesting of a significant mixing spreading offshore. However average 

346 values of P in Fig. [8] decrease when averaging over a wider area. 



347 The same inverse relationship is observed in Fig, 10 using chlorophyll data from Sea- 



348 nWIFS. This analysis confirms the result obtained from satellite velocity fields by Rossi 



D 

349 



et al. (2008, 2009) but using FSLEs computed on simulated velocity field with ROMS, at 



350 two different resolutions. In this case, the value of the slopes are: -3.5, -3.4 and -4.7 for 

351 Satellitel/4, ROMSl/4 and ROMSl/12, respectively. The fact that ROMS velocity data 

352 do not necessarily match the dates of SeaWIFS may explain the larger discrepancy in the 



353 values of the correlation coefficient showed in Fig, 10 



354 Then, let us present a brief description of the seasonal behavior of the system. In 

355 Fig. [TT] we display the temporal evolution over one year of the spatial averages of FSLEs 

356 (upper plot) and P (bottom). A climatological average for the case of Satellitel/4 using 

357 three years of data. We observe that the seasonal increase in mixing activity (from May 

358 to September, roughly winter) is associated to a decrease of the simulated phytoplankton. 

359 This also illustrates the seasonal inhibiting effect that the mixing activity has on the 

360 phytoplankton dynamics in winter. Note that the seasonal variation of light is not taken 

361 into account in our model. However, the temporal variability of plankton in the Benguela 

362 is mainly driven by the varying activity of the coastal upwelling cells, reproduced by the 

363 function S. 

364 Finally, a few sensitivity analysis were done to clarify the role of the 2D advection and 

365 the biological reactions in the simulated plankton fields. For this, we performed virtual 

366 experiments to determine the effect of both processes taken separately. A simulation with 

367 only advection of a passive tracer (without any upwelling parameterization) is compared 

368 to a similar simulation adding the biological reaction terms. The ad vect ion-only case 

369 reproduces well the smaller tracer concentrations in the southern domain, whereas the 



370 advection-reaction case presents a more constant latitudinal profile (see Fig. 12). This 
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371 confirms that the main influence on the spatial distribution of phytoplankton in the south 

372 is 2D advection, with the biological dynamics playing a minor role. 

373 Also, the per capita growth rate of N over time (i.e. N'^dN/dt) was computed and 

374 averaged over the coastal area in each subsystem to test the mechanisms proposed by 



375 [Gruber et al] ( |201lD (see Fig. [13] for the ROMS 1/1 2 simulation). We found that the 

376 mean value for each subsystem, North and South, are —3 ■ 10^^ and —1 ■ 10^^ day~^, 

377 respectively. This confirms that nutrients are being lost toward the open ocean by simple 

378 2D advection almost four times more in the south than in the north. It has to be compared 

379 with the mixing activity being about three times higher in the south than in the north 

380 (Fig. [3]). The same behavior is also observed in the other two cases ROMS 1/4 and 

381 Satellite (not shown). Note also that the loss of nutrient appear to be maximal in the 

382 winter months (maximum mixing), although there is a slight decay in between the two 

383 subsystems. 

384 5. Conclusions 

385 This study is based on numerical analysis from a simple biological NPZ model coupled 

386 with different velocity fields (satellite and model) over the Benguela area. Although in a 

387 simple framework, a reduction of phytoplankton concentrations in the coastal upwelling 

388 for increasing mesoscale activity has been successfully simulated. Horizontal stirring was 

389 estimated by computing the FSLEs and was correlated negatively with chlorophyll stocks. 

390 Similar results are found, though not presented in this manuscript, for the primary pro- 

n 

391 duction, defined as the first term in Fp (Eq,5), i.e. PP = (3 P. Some recent 

\J kn + N 

392 observational and modelling studies proposed the "nutrient leakage" as a mechanism to 

393 explain this negative correlation. Here we argue that Lagrangian Coherent Structures, 

394 mainly mesoscale eddies and filaments, transport a significant fraction of the recently 

395 upwelled nutrients nearshore toward the open ocean before being efficiently used by the 

396 pelagic food web. Although some studies dealt with 3D effect, we have shown that 2D 
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397 advection processes seems to play an important role in this suppressive effect. Our anal- 

398 ysis suggest that the inhibiting effect of the mesoscale activity on the plankton occurs 

399 when the mixing reach high levels, as in the south Benguela. However, this effect is not 

400 dominant under certain levels of turbulence. We have also shown that the inhibiting effect 

401 of intense mixing is maximal during the winter months. It might indicate that planktonic 

402 ecosystems in oceanic regions with vigorous mesoscale dynamics can be, as a first ap- 

403 proximation, easily modeled just by including a realistic flow fleld. The small residence 

404 times of waters in the productive area will smooth out all the other neglected biological 

405 factors in interaction. However, these factors are required when modelling an oceanic 

406 regions with low mixing, associated with high residence time leading to the predominance 

407 of complex combinations of factors. 

408 Our flndings conflrm the unexpected role that mesoscale activity has on biogeochemical 

409 dynamics in the productive coastal upwelling. Strong vertical velocities arc known to be 

410 associated with these physical structures and they might have another direct effect by 

411 transporting downward rich nutrient waters below the euphotic zone. Further studies are 

412 needed such as 3D realistic modelling that take into account the strong vertical dynamics 

413 in upwelling regions to test the complete mechanisms involved. 
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Figure 2: Spatial distribution of time average of weekly FSLE maps in the Benguela region, a) Three 
years average using data set SateUitel/4; b) one year average using ROMS 1/4; c) one year average using 
ROMSl/12. The units of the colorbar are I /days. 
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Figure 3: Zonal average over coastal bands of the FSLE time averages from Fig. [2] as a function of 
latitude. Top) From the coast to 3 degrees offshore; bottom) to 6 degrees offshore. 
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Figure 4: Spatial distribution of annual EKE in the Benguela region, a) using velocity data from Satellite 
at spatial resolution 1/4° {SateUitel/4) b) using velocity data from ROMS at spatial resolution 1/4° 
{ROMSl/4) c) using velocity data from ROMS at spatial resolution 1/12° {ROMSl/12). The units of 
the colorbar are {cm/s)^ 
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a) Satellitel/4, April 20, 2002 




d) ROMS 1/4, April 20 



g) ROMSl/12, April 20 
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b) Satellitel/4, June 14, 2002 



e) ROMS 1/4, June 14 



h) ROMSl/12, June 14 
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c) Satellitel/4, December 6, 2002 f) ROMSl/4, December 6 



i) ROMSl/12, December 6 
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Figure 5: Snapshots of spatial distribution of phytoplankton concentration from the simulations: Left 
column) corresponding to the simulation using Satellitel/4 ; Middle column) ROMS 1/4; Right column) 
from ROMSl/12. Logarithmic scale is used to improve the visualization of the structures. The units for 
the colorbar are mg/nv' 
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Figure 6: Snapshots every 8 days of large (top 30%) values of FSLE superimposed on P concentrations 
calculated from ROMSl/12 in mg/m^. Logarithmic scale for phytoplankton concentrations is used to 
improve the visualization of the structures 
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Figure 7: Spatial distribution of the time average of phytoplankton concentrations: a) Three years average 
using Satellitel/4, b) One year average from ROMSl/4, c) One year average from ROMSl/12, d) Three 
years average of monthly SeaWIFS data. The units of the colorbar are rag/nr' . 



29 




Latitude -0- (degrees) Latitude -9- (degrees) 



Figure 8: Zonal mean, over a 3 degrees (left) and 6 degrees (right) width coastal band, of the time 
averages of modelled phytoplankton (upper plots) and derived from satellite (lower plots) plotted as a 
function of latitude. 



30 



0.4 




0.08 0.1 0.12 0.14 0.16 0.18 
Mean FSLE (1/day) 



Figure 9: Weekly values of spatial averages of phytoplankton versus weekly values of spatial averages 
of FSLE, where the average are over the North and South subareas of Benguela. a) SateUitel/4, b) 
ROMS 1/4 and c) ROMS 1/12 
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Figure 10: Monthly values of spatial averages of Chlorophyll from SeaWIFS data versus spatial average 
of FSLE, where the average are over the North and South subareas of Benguela. FSLE values are from 
a) SatelUtel/l b) ROMSl/4 and c) ROMSl/12. 
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Figure 11: Temporal evolution of horizontal mixing (Spatial average of FSLEs) for the three velocity 
data sets (top) . Temporal evolution of spatial averages of simulated phytoplankton for the three velocity 
data sets (bottom). 
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Figure 12: Normalised comparison of the time averages of a passive scalar (advection only) and of P 
(advection-reaction), as a function of latitude. 
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Figure 13: Time evolution of the spatial average of the per capita growth rate of nutrients for the 
ROMSl/12 case. 
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